###### CANADA PANEL ANALYSIS

# setwd()

###### FUNCTIONS, PACKAGES, PARAMETERS ###### 

library(foreign)
library(readstata13)
library(reshape2)
library(plm)
library(plyr)
library(dplyr)
library(arm)
library(lmtest)
library(ggplot2)

# function to assist in plots

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  require(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- plyr::ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  # Rename the "mean" column    
  datac <- plyr::rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}



# default aesthetics for plots
pointsize<-4
linesize<-1
bartextsize<-4


cpanel <- read.csv("CanadaReplicationData.csv")


#####  DEFINE AND ELIMINATE THE WAVE-SPECIFIC RESPONDENTS
wave2.newrespondents <- wave2$responseid[which(!wave2$responseid%in%wave1$responseid)]
wave3.newrespondents <- wave3$responseid[which(!wave3$responseid%in%wave1$responseid)]
wave4.newrespondents <- wave4$responseid[which(!wave4$responseid%in%wave1$responseid)]
# no wave-specific respondents in wave 5

wave2.pool <- cpanel[which(cpanel$responseid%in%wave2.newrespondents),]
wave3.pool <- cpanel[which(cpanel$responseid%in%wave3.newrespondents),]
wave4.pool <- cpanel[which(cpanel$responseid%in%wave4.newrespondents),]

cpanel <- cpanel[-which(cpanel$responseid%in%wave2.newrespondents),]  
cpanel <- cpanel[-which(cpanel$responseid%in%wave3.newrespondents),]  
cpanel <- cpanel[-which(cpanel$responseid%in%wave4.newrespondents),]  



##### CREATE INDICATORS FOR COMPLETE WAVES, AND SUBSET

wave1_IDs <- unique(cpanel$responseid[which(cpanel$wave=="wave1")])
wave2_IDs <- unique(cpanel$responseid[which(cpanel$wave=="wave2")])
wave3_IDs <- unique(cpanel$responseid[which(cpanel$wave=="wave3")])
wave4_IDs <- unique(cpanel$responseid[which(cpanel$wave=="wave4")])
wave5_IDs <- unique(cpanel$responseid[which(cpanel$wave=="wave5")])

cpanel.complete5 <- cpanel[which(cpanel$responseid%in%wave5_IDs),]
cpanel.complete4 <- cpanel[which(cpanel$responseid%in%wave4_IDs),]
cpanel.complete3 <- cpanel[which(cpanel$responseid%in%wave3_IDs),]
cpanel.complete2 <- cpanel[which(cpanel$responseid%in%wave2_IDs),]






##### FIGURE 1: CARBON PRICING SUPPORT OVER TIME BY PROVINCE, 5 WAVES  #####

support_summary <- cpanel.complete5 %>%
  group_by(monthindex) %>%
  dplyr::summarize(supportBC = mean(support_cp[which(prov_clean=="BC")], na.rm=TRUE),
                   supportAB = mean(support_cp[which(prov_clean=="AB")], na.rm=TRUE),
                   supportSK = mean(support_cp[which(prov_clean=="SK")], na.rm=TRUE),
                   supportON = mean(support_cp[which(prov_clean=="ON")], na.rm=TRUE),
                   supportQC = mean(support_cp[which(prov_clean=="QC")], na.rm=TRUE),
                   seBC = sqrt(var(support_cp[which(prov_clean=="BC")], na.rm=TRUE)/length(support_cp[which(prov_clean=="BC")])),
                   seAB = sqrt(var(support_cp[which(prov_clean=="AB")], na.rm=TRUE)/length(support_cp[which(prov_clean=="AB")])),
                   seSK = sqrt(var(support_cp[which(prov_clean=="SK")], na.rm=TRUE)/length(support_cp[which(prov_clean=="SK")])),
                   seON = sqrt(var(support_cp[which(prov_clean=="ON")], na.rm=TRUE)/length(support_cp[which(prov_clean=="ON")])),
                   seQC = sqrt(var(support_cp[which(prov_clean=="QC")], na.rm=TRUE)/length(support_cp[which(prov_clean=="QC")])),
                   opposeBC = mean(oppose_cp[which(prov_clean=="BC")], na.rm=TRUE),
                   opposeAB = mean(oppose_cp[which(prov_clean=="AB")], na.rm=TRUE),
                   opposeSK = mean(oppose_cp[which(prov_clean=="SK")], na.rm=TRUE),
                   opposeON = mean(oppose_cp[which(prov_clean=="ON")], na.rm=TRUE),
                   opposeQC = mean(oppose_cp[which(prov_clean=="QC")], na.rm=TRUE))

supportmelt <- melt(support_summary, id="monthindex")
supportmelt$monthindex<- as.numeric(as.character(supportmelt$monthindex))

supportonly <- supportmelt[1:25,]
supportonly$se <- supportmelt$value[26:50]
quartz(width=7, height=5.5)

ggplot(supportonly, aes(x = monthindex, y = value*100)) + 
  geom_line(aes(color = variable,linetype=variable), size = linesize) + # aesthetics of the line between points
  geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
  scale_linetype_manual(name="",
                        values=c("solid", "solid", "solid", "solid", "solid"),# "dotted", "dotted", "dotted", "dotted", "dotted"),
                        labels=c("BC", "AB", "SK", "ON", "QC"))+
  scale_color_manual(name="",
                     values=c("#E69F00", "#56B4E9", "#009E73", "#D55E00","#CC79A7"),#, , "#D55E00", "#CC79A7")
                     labels=c("BC", "AB", "SK", "ON", "QC"))+#,
  geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
  scale_y_continuous(name="",expand = c(0,0),
                     limits = c(0,75),
                     breaks=c(0,25,50,75), labels=c("0%", "25%","50%", "75%")) + # controls the y axis
  scale_x_continuous(name="",breaks=c(1,3,6,10,14), labels=c("W1 (Feb)", "W2 (Apr)", "W3 (Jul)", "W4 (Nov)", "W5 (Apr)")) + # controls the x axis, including year ticks
  geom_vline(xintercept = 1.75, linetype="dotted")+
  geom_vline(xintercept = 4, linetype="solid")+
  geom_vline(xintercept = 9.5, linetype="dashed")+
  theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
        legend.text = element_text(size=9), # size of legend text
        legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
        legend.key=element_blank(), 
        aspect.ratio = 0.75, # of the entire grid (y vs. x axes)
        axis.text=element_text(size=9), # of the year ticks on the x axis
        axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
  theme( # remove the vertical grid lines
    panel.grid=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"))+ # add x and y axis lines
  geom_hline(yintercept=50,colour="gray")

 dev.copy2pdf(file="PaperFiguresFinal/Figure1.pdf")


 
 ##### FIGURE 2: EVALUATE SUPPORT BY PARTY IN REBATE VS. NON REBATE PROVINCES #####
 
 support_summary <- cpanel.complete5 %>%
   group_by(monthindex) %>%
   dplyr::summarize(supportLIBnorebate = mean(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="BC"|prov_clean=="QC"))], na.rm=TRUE),
                    supportLIBrebate = mean(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="SK"|prov_clean=="ON"))], na.rm=TRUE),
                    supportCONnorebate = mean(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="BC"|prov_clean=="QC"))], na.rm=TRUE),
                    supportCONrebate = mean(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="SK"|prov_clean=="ON"))], na.rm=TRUE),
                    seLIBnorebate = sqrt(var(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="BC"|prov_clean=="QC"))], 
                                             na.rm=TRUE)/length(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="BC"|prov_clean=="QC"))])),
                    seLIBrebate = sqrt(var(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="SK"|prov_clean=="ON"))], 
                                           na.rm=TRUE)/length(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="SK"|prov_clean=="ON"))])),
                    seCONnorebate = sqrt(var(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="BC"|prov_clean=="QC"))], 
                                             na.rm=TRUE)/length(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="BC"|prov_clean=="QC"))])),
                    seCONrebate = sqrt(var(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="SK"|prov_clean=="ON"))], 
                                           na.rm=TRUE)/length(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="SK"|prov_clean=="ON"))])))
 
 supportmelt <- melt(support_summary, id="monthindex")
 supportmelt$monthindex<- as.numeric(as.character(supportmelt$monthindex))
 
 supportonly <- supportmelt[1:20,]
 supportonly$se <- supportmelt$value[21:40]
 quartz(width=7, height=5.5)
 
 ggplot(supportonly, aes(x = monthindex, y = value*100)) + 
   geom_line(aes(color = variable,linetype=variable), size = linesize) + # aesthetics of the line between points
   geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
   scale_linetype_manual(name="",
                         values=c("solid", "solid", "solid", "solid"),# "dotted", "dotted", "dotted", "dotted", "dotted"),
                         labels=c("Liberal (QC, BC)", "Liberal (SK, ON)", "Conservative (QC, BC)", "Conservative (SK, ON)"))+
   scale_color_manual(name="",
                      values=c("#E69F00", "#56B4E9", "#009E73", "#D55E00"),#"#F0E442"),#, , "#D55E00", "#CC79A7")
                      labels=c("Liberal (QC, BC)", "Liberal (SK, ON)", "Conservative (QC, BC)", "Conservative (SK, ON)"))+#,
   geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
   scale_y_continuous(name="",expand = c(0,0),
                      limits = c(0,100),
                      breaks=c(0,25,50,75,100), labels=c("0%", "25%","50%", "75%", "100%")) + # controls the y axis
   scale_x_continuous(name="",breaks=c(1,3,6,10,14), labels=c("W1 (Feb)", "W2 (Apr)", "W3 (Jul)", "W4 (Nov)", "W5 (Apr)")) + # controls the x axis, including year ticks
   geom_vline(xintercept = 1.75, linetype="dotted")+
   geom_vline(xintercept = 4, linetype="solid")+
   geom_vline(xintercept = 9.5, linetype="dashed")+
   theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
         legend.position="bottom", legend.direction = "vertical",# moves legend to the top and makes it stacked
         legend.text = element_text(size=9), # size of legend text
         legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
         legend.key=element_blank(), 
         aspect.ratio = 0.75, # of the entire grid (y vs. x axes)
         axis.text=element_text(size=9), # of the year ticks on the x axis
         axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
   
   theme( # remove the vertical grid lines
     panel.grid=element_blank(),
    
     panel.background = element_blank(),
     panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
     axis.line=element_line(colour="gray"))+ # add x and y axis lines
   geom_hline(yintercept=50,colour="gray")
 
 dev.copy2pdf(file="PaperFiguresFinal/Figure2.pdf")
 
 
 ##### TABLE 1 AND ACCOMPANYING PARAGRAPHS #####
 
# Table 1 Data (Perceived Rebate Amount by Province)
 mean(cpanel$divestimate[which(cpanel$prov_clean=="BC"&cpanel$wave=="wave3")], na.rm=TRUE)
 mean(cpanel$divestimate[which(cpanel$prov_clean=="AB"&cpanel$wave=="wave3")], na.rm=TRUE)
 mean(cpanel$divestimate[which(cpanel$prov_clean=="QC"&cpanel$wave=="wave3")], na.rm=TRUE)
 sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="BC"&cpanel$wave=="wave3")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="BC"&cpanel$wave=="wave3")]))
 sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="AB"&cpanel$wave=="wave3")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="AB"&cpanel$wave=="wave3")]))
 sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="QC"&cpanel$wave=="wave3")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="QC"&cpanel$wave=="wave3")]))
 
 mean(cpanel$divestimate[which(cpanel$prov_clean=="ON"&cpanel$wave=="wave3")], na.rm=TRUE)
 mean(cpanel$divestimate[which(cpanel$prov_clean=="SK"&cpanel$wave=="wave3")], na.rm=TRUE)
 sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="ON"&cpanel$wave=="wave3")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="ON"&cpanel$wave=="wave3")]))
 sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="SK"&cpanel$wave=="wave3")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="SK"&cpanel$wave=="wave3")]))
 
 # Table 1 Data (True Rebate Amounts by Province  )
 mean(cpanel$federal_realdividend_raw[which(cpanel$prov_clean=="ON"&cpanel$wave=="wave3")], na.rm=TRUE)
 mean(cpanel$federal_realdividend_raw[which(cpanel$prov_clean=="SK"&cpanel$wave=="wave3")], na.rm=TRUE)
 
 
# Perceived rebate amounts subset to individuals who *think* they received a rebate
 mean(cpanel$divestimate[which(cpanel$prov_clean=="ON"&cpanel$think_rebate=="Yes"&cpanel$wave=="wave3")], na.rm=TRUE)
 mean(cpanel$divestimate[which(cpanel$prov_clean=="SK"&cpanel$think_rebate=="Yes"&cpanel$wave=="wave3")], na.rm=TRUE)
 sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="ON"&cpanel$think_rebate=="Yes"&cpanel$wave=="wave3")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="ON"&cpanel$think_rebate=="Yes"&cpanel$wave=="wave3")]))
 sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="SK"&cpanel$think_rebate=="Yes"&cpanel$wave=="wave3")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="SK"&cpanel$think_rebate=="Yes"&cpanel$wave=="wave3")]))
 

##### FIGURE 3: SEE REPLICATION CODE FOR SWISS SURVEY #####
 

##### FIGURE 4: EXPERIMENTAL RESULTS AND COEFFICIENT PLOT #####

#subset to experimental cases
exp_analysis <- cpanel[which(cpanel$wave=="wave4"),]
exp_analysis<- exp_analysis[which(exp_analysis$prov_clean=="ON"|exp_analysis$prov_clean=="SK"),] 


# results reported in text
exp_results_support <- lm(supportDV~treat, data=exp_analysis)
summary(exp_results_support)

exp_results_netcost <- lm(netcostdv~treat, data=exp_analysis)
summary(exp_results_netcost)

exp_results_netcost <- lm(netcostdv~treat, data=exp_analysis[which(exp_analysis$wave4_party=="Conservatives"),])
summary(exp_results_netcost)


# reproduce Figure  4 (top)
panel_cpsupport <- summarySE(data=exp_analysis, measurevar=c("supportDV"), groupvars=c("treat"), na.rm=TRUE)
panel_cpsupportLIB <- summarySE(data=exp_analysis[which(exp_analysis$wave4_party=="Liberals"),], measurevar=c("supportDV"), groupvars=c("treat"), na.rm=TRUE)
panel_cpsupportNDP <- summarySE(data=exp_analysis[which(exp_analysis$wave4_party=="ndp"),], measurevar=c("supportDV"), groupvars=c("treat"), na.rm=TRUE)
panel_cpsupportCON <- summarySE(data=exp_analysis[which(exp_analysis$wave4_party=="Conservatives"),], measurevar=c("supportDV"), groupvars=c("treat"), na.rm=TRUE)

canada_cpsupportgraph <- rbind (panel_cpsupport, panel_cpsupportLIB, panel_cpsupportNDP, panel_cpsupportCON)
canada_cpsupportgraph$group <- c("all", "all", "Liberal", "Liberal", "NDP", "NDP", "Conservative", "Conservative")

quartz(width=7, height=5.5)

ggplot(canada_cpsupportgraph, aes(x=as.factor(treat), y=supportDV, color=group)) + 
  geom_point(position=position_dodge(.4), size=4)+
  scale_y_continuous(limits=c(1,4)) +
  scale_color_manual(values=c("#000000","#56B4E9","#009E73","#D55E00"))+
  ylab("Support for existing policy (4 point scale)") +
  xlab("Treatment Status") +
  theme_bw()+
  geom_errorbar(aes(ymin=supportDV-ci, ymax=supportDV+ci),
                width=.2,
                size=1,
                position=position_dodge(.4))


dev.copy2pdf(file="PaperFiguresFinal/Figure4_top.pdf")


# reproduce Figure  4 (bottom)

panel_cpcost <- summarySE(data=exp_analysis, measurevar=c("netcostdv"), groupvars=c("treat"), na.rm=TRUE)
panel_cpcostLIB <- summarySE(data=exp_analysis[which(exp_analysis$wave4_party=="Liberals"),], measurevar=c("netcostdv"), groupvars=c("treat"), na.rm=TRUE)
panel_cpcostNDP <- summarySE(data=exp_analysis[which(exp_analysis$wave4_party=="ndp"),], measurevar=c("netcostdv"), groupvars=c("treat"), na.rm=TRUE)
panel_cpcostCON <- summarySE(data=exp_analysis[which(exp_analysis$wave4_party=="Conservatives"),], measurevar=c("netcostdv"), groupvars=c("treat"), na.rm=TRUE)

canada_cpcostgraph <- rbind (panel_cpcost, panel_cpcostLIB, panel_cpcostNDP, panel_cpcostCON)
canada_cpcostgraph$group <- c("all", "all", "Liberal", "Liberal", "NDP", "NDP", "Conservative", "Conservative")

quartz(width=7, height=5.5)

ggplot(canada_cpcostgraph, aes(x=as.factor(treat), y=netcostdv, color=group)) + 
  geom_point(position=position_dodge(.4), size=4)+
  scale_y_continuous(limits=c(-1,1)) +
  scale_color_manual(values=c("#000000","#56B4E9","#009E73","#D55E00"))+
  ylab("Belief receive more than pay (scale from -1 to 1)") +
  xlab("Treatment Status") +
  theme_bw()+
  geom_errorbar(aes(ymin=netcostdv-ci, ymax=netcostdv+ci),
                width=.2,       
                size=1,
                position=position_dodge(.4))

dev.copy2pdf(file="PaperFiguresFinal/Figure4_bottom.pdf")

##### FIGURE 5: SEE REPLICATION CODE FOR SWISS SURVEY #####


##### EXTENDED DATA FIGURE 1: CARBON PRICING OPPOSITION BY PROVINCE AND SURVEY WAVE #####

# produce any oppose over time figure for the 5 wave panel
oppose_summary <- cpanel.complete5 %>%
  group_by(monthindex) %>%
  dplyr::summarize(supportBC = mean(oppose_cp[which(prov_clean=="BC")], na.rm=TRUE),
                   supportAB = mean(oppose_cp[which(prov_clean=="AB")], na.rm=TRUE),
                   supportSK = mean(oppose_cp[which(prov_clean=="SK")], na.rm=TRUE),
                   supportON = mean(oppose_cp[which(prov_clean=="ON")], na.rm=TRUE),
                   supportQC = mean(oppose_cp[which(prov_clean=="QC")], na.rm=TRUE),
                   seBC = sqrt(var(oppose_cp[which(prov_clean=="BC")], na.rm=TRUE)/length(oppose_cp[which(prov_clean=="BC")])),
                   seAB = sqrt(var(oppose_cp[which(prov_clean=="AB")], na.rm=TRUE)/length(oppose_cp[which(prov_clean=="AB")])),
                   seSK = sqrt(var(oppose_cp[which(prov_clean=="SK")], na.rm=TRUE)/length(oppose_cp[which(prov_clean=="SK")])),
                   seON = sqrt(var(oppose_cp[which(prov_clean=="ON")], na.rm=TRUE)/length(oppose_cp[which(prov_clean=="ON")])),
                   seQC = sqrt(var(oppose_cp[which(prov_clean=="QC")], na.rm=TRUE)/length(oppose_cp[which(prov_clean=="QC")])))

opposemelt <- melt(oppose_summary, id="monthindex")
opposemelt$monthindex<- as.numeric(as.character(opposemelt$monthindex))

opposeonly <- opposemelt[1:25,]
opposeonly$se <- opposemelt$value[26:50]

quartz(width=7, height=5.5)

ggplot(opposeonly, aes(x = monthindex, y = value*100)) + 
  geom_line(aes(color = variable,linetype=variable), size = linesize) + # aesthetics of the line between points
  geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
  scale_linetype_manual(name="",
                        values=c("solid", "solid", "solid", "solid", "solid"),# "dotted", "dotted", "dotted", "dotted", "dotted"),
                        labels=c("BC", "AB", "SK", "ON", "QC"))+
  scale_color_manual(name="",
                     values=c("#E69F00", "#56B4E9", "#009E73", "#D55E00","#CC79A7"),#, , "#D55E00", "#CC79A7")
                     labels=c("BC", "AB", "SK", "ON", "QC"))+#,
  geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
  scale_y_continuous(name="",expand = c(0,0),
                     limits = c(0,75),
                     breaks=c(0,25,50,75), labels=c("0%", "25%","50%", "75%")) + # controls the y axis
  scale_x_continuous(name="",breaks=c(1,3,6,10,14), labels=c("W1 (Feb)", "W2 (Apr)", "W3 (Jul)", "W4 (Nov)", "W5 (Apr)")) + # controls the x axis, including year ticks
  geom_vline(xintercept = 1.75, linetype="dotted")+
  geom_vline(xintercept = 4, linetype="solid")+
  geom_vline(xintercept = 9.5, linetype="dashed")+
  theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
        legend.text = element_text(size=9), # size of legend text
        legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
        legend.key=element_blank(), 
        aspect.ratio = 0.75, # of the entire grid (y vs. x axes)
        axis.text=element_text(size=9), # of the year ticks on the x axis
        axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
  
  theme( # remove the vertical grid lines
    panel.grid=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"))+ # add x and y axis lines
  geom_hline(yintercept=50,colour="gray")

dev.copy2pdf(file="PaperFiguresFinal/ExtendedDataFigure1.pdf")


##### EXTENDED DATA FIGURE 2: CARBON PRICING SUPPORT BY REBATE AND PARTY PRESENCE (INCLUDING ALBERTA) #####


# also, produce robustness figure for the 5 wave panel and party preferences
support_summary <- cpanel.complete5 %>%
  group_by(monthindex) %>%
  dplyr::summarize(supportLIBnorebate = mean(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="BC"|prov_clean=="QC"|prov_clean=="AB"))], na.rm=TRUE),
                   supportLIBrebate = mean(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="SK"|prov_clean=="ON"))], na.rm=TRUE),
                   supportCONnorebate = mean(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="BC"|prov_clean=="QC"|prov_clean=="AB"))], na.rm=TRUE),
                   supportCONrebate = mean(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="SK"|prov_clean=="ON"))], na.rm=TRUE),
                   seLIBnorebate = sqrt(var(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="BC"|prov_clean=="QC"|prov_clean=="AB"))], 
                                            na.rm=TRUE)/length(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="BC"|prov_clean=="QC"|prov_clean=="AB"))])),
                   seLIBrebate = sqrt(var(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="SK"|prov_clean=="ON"))], 
                                          na.rm=TRUE)/length(support_cp[which(wave1_party=="Liberal Party"&(prov_clean=="SK"|prov_clean=="ON"))])),
                   seCONnorebate = sqrt(var(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="BC"|prov_clean=="QC"|prov_clean=="AB"))], 
                                            na.rm=TRUE)/length(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="BC"|prov_clean=="QC"|prov_clean=="AB"))])),
                   seCONrebate = sqrt(var(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="SK"|prov_clean=="ON"))], 
                                          na.rm=TRUE)/length(support_cp[which(wave1_party=="Conservative Party"&(prov_clean=="SK"|prov_clean=="ON"))])))

supportmelt <- melt(support_summary, id="monthindex")
supportmelt$monthindex<- as.numeric(as.character(supportmelt$monthindex))

supportonly <- supportmelt[1:20,]
supportonly$se <- supportmelt$value[21:40]
quartz(width=7, height=5.5)

ggplot(supportonly, aes(x = monthindex, y = value*100)) + 
  geom_line(aes(color = variable,linetype=variable), size = linesize) + # aesthetics of the line between points
  geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
  scale_linetype_manual(name="",
                        values=c("solid", "solid", "solid", "solid"),# "dotted", "dotted", "dotted", "dotted", "dotted"),
                        labels=c("Liberal (QC, BC, AB)", "Liberal (SK, ON)", "Conservative (QC, BC, AB)", "Conservative (SK, ON)"))+
  #  "Oppose (BC)", "Oppose (AB)", "Oppose (SK)", "Oppose (ON)", "Oppose (QC)"))+
  scale_color_manual(name="",
                     values=c("#E69F00", "#56B4E9", "#009E73", "#D55E00"),#"#F0E442"),#, , "#D55E00", "#CC79A7")
                     labels=c("Liberal (QC, BC, AB)", "Liberal (SK, ON)", "Conservative (QC, BC, AB)", "Conservative (SK, ON)"))+#,
  geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
  scale_y_continuous(name="",expand = c(0,0),
                     limits = c(0,100),
                     breaks=c(0,25,50,75,100), labels=c("0%", "25%","50%", "75%", "100%")) + # controls the y axis
  scale_x_continuous(name="",breaks=c(1,3,6,10,14), labels=c("W1 (Feb)", "W2 (Apr)", "W3 (Jul)", "W4 (Nov)", "W5 (Apr)")) + # controls the x axis, including year ticks
  geom_vline(xintercept = 1.75, linetype="dotted")+
  geom_vline(xintercept = 4, linetype="solid")+
  geom_vline(xintercept = 9.5, linetype="dashed")+
  theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
        legend.position="bottom", legend.direction = "vertical",# moves legend to the top and makes it stacked
        legend.text = element_text(size=9), # size of legend text
        legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
        legend.key=element_blank(), 
        aspect.ratio = 0.75, # of the entire grid (y vs. x axes)
        axis.text=element_text(size=9), # of the year ticks on the x axis
        axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
  
  theme( # remove the vertical grid lines
    panel.grid=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"))+ # add x and y axis lines
  geom_hline(yintercept=50,colour="gray")

dev.copy2pdf(file="PaperFiguresFinal/ExtendedDataFigure3.pdf")


##### EXTENDED DATA FIGURE 3: CARBON PRICING SUPPORT BY PARTY AND COST EXPOSURE #####

support_summary <- cpanel.complete5 %>%
  group_by(monthindex) %>%
  dplyr::summarize(supportLIBdrivealone = mean(support_cp[which(wave1_party=="Liberal Party"&(wave1_drive=="Drive alone"&wave1_drive!="This doesn't apply to me"))], na.rm=TRUE),
                   supportLIBother = mean(support_cp[which(wave1_party=="Liberal Party"&(wave1_drive!="Drive alone"&wave1_drive!="This doesn't apply to me"))], na.rm=TRUE),
                   supportCONdrivealone = mean(support_cp[which(wave1_party=="Conservative Party"&(wave1_drive=="Drive alone"&wave1_drive!="This doesn't apply to me"))], na.rm=TRUE),
                   supportCONother = mean(support_cp[which(wave1_party=="Conservative Party"&(wave1_drive!="Drive alone"&wave1_drive!="This doesn't apply to me"))], na.rm=TRUE),
                   seLIBcommute = sqrt(var(support_cp[which(wave1_party=="Liberal Party"&(wave1_drive=="Drive alone"&wave1_drive!="This doesn't apply to me"))], 
                                           na.rm=TRUE)/length(support_cp[which(wave1_party=="Liberal Party"&(wave1_drive=="Drive alone"&wave1_drive!="This doesn't apply to me"))])),
                   seLIBnocommute = sqrt(var(support_cp[which(wave1_party=="Liberal Party"&(wave1_drive!="Drive alone"&wave1_drive!="This doesn't apply to me"))], 
                                             na.rm=TRUE)/length(support_cp[which(wave1_party=="Liberal Party"&(wave1_drive!="Drive alone"&wave1_drive!="This doesn't apply to me"))])),
                   seCONcommute = sqrt(var(support_cp[which(wave1_party=="Conservative Party"&(wave1_drive=="Drive alone"&wave1_drive!="This doesn't apply to me"))], 
                                           na.rm=TRUE)/length(support_cp[which(wave1_party=="Conservative Party"&(wave1_drive=="Drive alone"&wave1_drive!="This doesn't apply to me"))])),
                   seCONnocommute = sqrt(var(support_cp[which(wave1_party=="Conservative Party"&(wave1_drive!="Drive alone"&wave1_drive!="This doesn't apply to me"))], 
                                             na.rm=TRUE)/length(support_cp[which(wave1_party=="Conservative Party"&(wave1_drive!="Drive alone"&wave1_drive!="This doesn't apply to me"))])))

supportmelt <- melt(support_summary, id="monthindex")
supportmelt$monthindex<- as.numeric(as.character(supportmelt$monthindex))

supportonly <- supportmelt[1:20,]
supportonly$se <- supportmelt$value[21:40]

quartz(width=7, height=5.5)

ggplot(supportonly, aes(x = monthindex, y = value*100)) + 
  geom_line(aes(color = variable,linetype=variable), size = linesize) + # aesthetics of the line between points
  geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
  scale_linetype_manual(name="",
                        values=c("solid", "solid", "solid", "solid"),# "dotted", "dotted", "dotted", "dotted", "dotted"),
                        labels=c("Liberals (drive alone)", "Liberals (other)", "Conservatives (drive alone)", "Conservatives (other)"))+#,
  scale_color_manual(name="",
                     values=c("#E69F00", "#56B4E9", "#009E73", "#D55E00"),#"#F0E442"),#, , "#D55E00", "#CC79A7")
                     labels=c("Liberals (drive alone)", "Liberals (other)", "Conservatives (drive alone)", "Conservatives (other)"))+#,
  geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
  scale_y_continuous(name="",expand = c(0,0),
                     limits = c(0,100),
                     breaks=c(0,25,50,75,100), labels=c("0%", "25%","50%", "75%", "100%")) + # controls the y axis
  scale_x_continuous(name="",breaks=c(1,3,6,10,14), labels=c("W1 (Feb)", "W2 (Apr)", "W3 (Jul)", "W4 (Nov)", "W5 (Apr)")) + # controls the x axis, including year ticks
  geom_vline(xintercept = 1.75, linetype="dotted")+
  geom_vline(xintercept = 4, linetype="solid")+
  geom_vline(xintercept = 9.5, linetype="dashed")+
  theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
        legend.position="bottom", legend.direction = "vertical",# moves legend to the top and makes it stacked
        legend.text = element_text(size=9), # size of legend text
        legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
        legend.key=element_blank(), 
        aspect.ratio = 0.75, # of the entire grid (y vs. x axes)
        axis.text=element_text(size=9), # of the year ticks on the x axis
        axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
  
  theme( # remove the vertical grid lines
    panel.grid=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"))+ # add x and y axis lines
  geom_hline(yintercept=50,colour="gray")

dev.copy2pdf(file="PaperFiguresFinal/ExtendedDataFigure3.pdf")



##### EXTENDED DATA FIGURE 4: DISTRIBUTION OF PERCEIVED REBATE SIZES IN CANADA #####

divcheck <- cpanel[which(cpanel$prov_clean=="ON"|cpanel$prov_clean=="SK"),]
divcheck <- divcheck[which(divcheck$wave=="wave3"),]

real100 <- divcheck[which(divcheck$federal_realdividend_raw>100&divcheck$federal_realdividend_raw<200),] #"$100-200"
real200 <- divcheck[which(divcheck$federal_realdividend_raw>200&divcheck$federal_realdividend_raw<300),] #"$200-300"
real300 <- divcheck[which(divcheck$federal_realdividend_raw>300&divcheck$federal_realdividend_raw<400),] #"$300-400"
real400 <- divcheck[which(divcheck$federal_realdividend_raw>400&divcheck$federal_realdividend_raw<500),] #"$400-500"
real500 <- divcheck[which(divcheck$federal_realdividend_raw>500&divcheck$federal_realdividend_raw<600),] #"$500-600"
real600 <- divcheck[which(divcheck$federal_realdividend_raw>600&divcheck$federal_realdividend_raw<700),] #"$600-700"
real700 <- divcheck[which(divcheck$federal_realdividend_raw>700&divcheck$federal_realdividend_raw<800),] #"$700-800"

real100$claimd <- factor(real100$claimd, levels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"))
real200$claimd <- factor(real200$claimd, levels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"))
real300$claimd <- factor(real300$claimd, levels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"))
real400$claimd <- factor(real400$claimd, levels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"))
real500$claimd <- factor(real500$claimd, levels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"))
real600$claimd <- factor(real600$claimd, levels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"))
real700$claimd <- factor(real700$claimd, levels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"))


ggplot(real100, aes(x=claimd, fill=claimd)) + 
  geom_bar(show.legend=FALSE)+
  scale_y_continuous(limits=c(0,50)) +
  scale_x_discrete(labels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"), drop=FALSE)+
  scale_fill_manual(values=c("#000000","#000000","#009E73", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000"))+
  ylab("Number of survey respondents") +
  xlab("Perceived household rebate range") +
  ggtitle("Actual received rebate between $100 and $200")+
  theme_bw()
dev.copy2pdf(file="PaperFigures/ExtendedDataFigure4_topleft.pdf")

ggplot(real200, aes(x=claimd, fill=claimd)) + 
  geom_bar(show.legend=FALSE)+
  scale_y_continuous(limits=c(0,50)) +
  scale_x_discrete(labels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"), drop=FALSE)+
  scale_fill_manual(values=c("#000000","#000000","#000000", "#009E73", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000"))+
  ylab("Number of survey respondents") +
  xlab("Perceived household rebate range") +
  ggtitle("Actual received rebate between $200 and $300")+
  theme_bw()
dev.copy2pdf(file="PaperFigures/ExtendedDataFigure4_topright.pdf")

ggplot(real300, aes(x=claimd, fill=claimd)) + 
  geom_bar(show.legend=FALSE)+
  scale_y_continuous(limits=c(0,50)) +
  scale_x_discrete(labels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"), drop=FALSE)+
  scale_fill_manual(values=c("#000000","#000000","#000000", "#000000", "#009E73", "#000000", "#000000", "#000000", "#000000", "#000000", "#000000"), drop=FALSE)+
  ylab("Number of survey respondents") +
  xlab("Perceived household rebate range") +
  ggtitle("Actual received rebate between $300 and $400")+
  theme_bw()
dev.copy2pdf(file="PaperFigures/ExtendedDataFigure4_centerleft.pdf")

ggplot(real400, aes(x=claimd, fill=claimd)) + 
  geom_bar(show.legend=FALSE)+
  scale_y_continuous(limits=c(0,50)) +
  scale_x_discrete(labels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"), drop=FALSE)+
  scale_fill_manual(values=c("#000000","#000000","#000000", "#000000", "#000000", "#009E73", "#000000", "#000000", "#000000", "#000000", "#000000"), drop=FALSE)+
  ylab("Number of survey respondents") +
  xlab("Perceived household rebate range") +
  ggtitle("Actual received rebate between $400 and $500")+
  theme_bw()
dev.copy2pdf(file="PaperFigures/ExtendedDataFigure4_centerright.pdf")

rebate500 <- ggplot(real500, aes(x=claimd, fill=claimd)) + 
  geom_bar(show.legend=FALSE)+
  scale_y_continuous(limits=c(0,50)) +
  scale_x_discrete(labels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"), drop=FALSE)+
  scale_fill_manual(values=c("#000000","#000000","#000000", "#000000", "#000000", "#000000", "#009E73", "#000000", "#000000", "#000000", "#000000"), drop=FALSE)+
  ylab("Number of survey respondents") +
  xlab("Perceived household rebate range") +
  ggtitle("Actual received rebate between $500 and $600")+
  theme_bw()
dev.copy2pdf(file="PaperFigures/ExtendedDataFigure4_bottomleft.pdf")

rebate600 <- ggplot(real600, aes(x=claimd, fill=claimd)) + 
  geom_bar(show.legend=FALSE)+
  scale_y_continuous(limits=c(0,50)) +
  scale_x_discrete(labels=c("$0", "$1-100", "$100-200", "$200-300","$300-400","$400-500","$500-600","$600-700","$700-800", "$800-900","$1000+"), drop=FALSE)+
  scale_fill_manual(values=c("#000000","#000000","#000000", "#000000", "#000000", "#000000", "#009E73", "#000000", "#000000", "#000000", "#000000"), drop=FALSE)+
  ylab("Number of survey respondents") +
  xlab("Perceived household rebate range") +
  ggtitle("Actual received rebate between $600 and $700")+
  theme_bw()
dev.copy2pdf(file="PaperFigures/ExtendedDataFigure4_bottomright.pdf")


##### EXTENDED DATA FIGURES 5,6,7,8: No code
 
##### SI SECTION 1: No code ##### 

##### SI SECTION 2: No code ##### 

#####  SI SECTION 3: CARBON PRICING SUPPORT ACROSS FIRST FOUR PANEL WAVES ##### 
support_summary <- cpanel.complete4 %>%
  group_by(monthindex) %>%
  dplyr::summarize(supportBC = mean(support_cp[which(prov_clean=="BC")], na.rm=TRUE),
                   supportAB = mean(support_cp[which(prov_clean=="AB")], na.rm=TRUE),
                   supportSK = mean(support_cp[which(prov_clean=="SK")], na.rm=TRUE),
                   supportON = mean(support_cp[which(prov_clean=="ON")], na.rm=TRUE),
                   supportQC = mean(support_cp[which(prov_clean=="QC")], na.rm=TRUE),
                   seBC = sqrt(var(support_cp[which(prov_clean=="BC")], na.rm=TRUE)/length(support_cp[which(prov_clean=="BC")])),
                   seAB = sqrt(var(support_cp[which(prov_clean=="AB")], na.rm=TRUE)/length(support_cp[which(prov_clean=="AB")])),
                   seSK = sqrt(var(support_cp[which(prov_clean=="SK")], na.rm=TRUE)/length(support_cp[which(prov_clean=="SK")])),
                   seON = sqrt(var(support_cp[which(prov_clean=="ON")], na.rm=TRUE)/length(support_cp[which(prov_clean=="ON")])),
                   seQC = sqrt(var(support_cp[which(prov_clean=="QC")], na.rm=TRUE)/length(support_cp[which(prov_clean=="QC")])),
                   opposeBC = mean(oppose_cp[which(prov_clean=="BC")], na.rm=TRUE),
                   opposeAB = mean(oppose_cp[which(prov_clean=="AB")], na.rm=TRUE),
                   opposeSK = mean(oppose_cp[which(prov_clean=="SK")], na.rm=TRUE),
                   opposeON = mean(oppose_cp[which(prov_clean=="ON")], na.rm=TRUE),
                   opposeQC = mean(oppose_cp[which(prov_clean=="QC")], na.rm=TRUE))

supportmelt <- melt(support_summary, id="monthindex")
supportmelt$monthindex<- as.numeric(as.character(supportmelt$monthindex))
supportmelt <- supportmelt [-which(supportmelt$monthindex==14),] # restrict to four waves for graph

supportonly <- supportmelt[1:20,]
supportonly$se <- supportmelt$value[21:40]

quartz(width=7, height=5.5)

ggplot(supportonly, aes(x = monthindex, y = value*100)) + 
  geom_line(aes(color = variable,linetype=variable), size = linesize) + # aesthetics of the line between points
  geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
  scale_linetype_manual(name="",
                        values=c("solid", "solid", "solid", "solid", "solid"),# "dotted", "dotted", "dotted", "dotted", "dotted"),
                        labels=c("BC", "AB", "SK", "ON", "QC"))+
  scale_color_manual(name="",
                     values=c("#E69F00", "#56B4E9", "#009E73", "#D55E00","#CC79A7"),#, , "#D55E00", "#CC79A7")
                     labels=c("BC", "AB", "SK", "ON", "QC"))+#,
  geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
  scale_y_continuous(name="",expand = c(0,0),
                     limits = c(0,75),
                     breaks=c(0,25,50,75), labels=c("0%", "25%","50%", "75%")) + # controls the y axis
  scale_x_continuous(name="",breaks=c(1,3,6,10), labels=c("W1 (Feb)", "W2 (Apr)", "W3 (Jul)", "W4 (Nov)")) + # controls the x axis, including year ticks
  geom_vline(xintercept = 1.75, linetype="dotted")+
  geom_vline(xintercept = 4, linetype="solid")+
  geom_vline(xintercept = 9.5, linetype="dashed")+
  theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
        legend.text = element_text(size=9), # size of legend text
        legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
        legend.key=element_blank(), 
        axis.text=element_text(size=9), # of the year ticks on the x axis
        axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
  
  theme( # remove the vertical grid lines
    panel.grid=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"))+ # add x and y axis lines
  geom_hline(yintercept=50,colour="gray")

dev.copy2pdf(file="PaperFiguresFinal/support_overtime4.pdf")



##### SI SECTION 4: STRONG CARBON PRICING OPPOSITION OVER TIME ##### 

# produce strong oppose over time figure for the 5 wave panel
strongoppose_summary <- cpanel.complete5 %>%
  group_by(monthindex) %>%
  dplyr::summarize(supportBC = mean(strongoppose_cp[which(prov_clean=="BC")], na.rm=TRUE),
                   supportAB = mean(strongoppose_cp[which(prov_clean=="AB")], na.rm=TRUE),
                   supportSK = mean(strongoppose_cp[which(prov_clean=="SK")], na.rm=TRUE),
                   supportON = mean(strongoppose_cp[which(prov_clean=="ON")], na.rm=TRUE),
                   supportQC = mean(strongoppose_cp[which(prov_clean=="QC")], na.rm=TRUE),
                   seBC = sqrt(var(strongoppose_cp[which(prov_clean=="BC")], na.rm=TRUE)/length(strongoppose_cp[which(prov_clean=="BC")])),
                   seAB = sqrt(var(strongoppose_cp[which(prov_clean=="AB")], na.rm=TRUE)/length(strongoppose_cp[which(prov_clean=="AB")])),
                   seSK = sqrt(var(strongoppose_cp[which(prov_clean=="SK")], na.rm=TRUE)/length(strongoppose_cp[which(prov_clean=="SK")])),
                   seON = sqrt(var(strongoppose_cp[which(prov_clean=="ON")], na.rm=TRUE)/length(strongoppose_cp[which(prov_clean=="ON")])),
                   seQC = sqrt(var(strongoppose_cp[which(prov_clean=="QC")], na.rm=TRUE)/length(strongoppose_cp[which(prov_clean=="QC")])))

strongopposemelt <- melt(strongoppose_summary, id="monthindex")
strongopposemelt$monthindex<- as.numeric(as.character(strongopposemelt$monthindex))

strongopposeonly <- strongopposemelt[1:25,]
strongopposeonly$se <- strongopposemelt$value[26:50]

quartz(width=7, height=5.5)

ggplot(strongopposeonly, aes(x = monthindex, y = value*100)) + 
  geom_line(aes(color = variable,linetype=variable), size = linesize) + # aesthetics of the line between points
  geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
  scale_linetype_manual(name="",
                        values=c("solid", "solid", "solid", "solid", "solid"),# "dotted", "dotted", "dotted", "dotted", "dotted"),
                        labels=c("BC", "AB", "SK", "ON", "QC"))+
  scale_color_manual(name="",
                     values=c("#E69F00", "#56B4E9", "#009E73", "#D55E00","#CC79A7"),#, , "#D55E00", "#CC79A7")
                     labels=c("BC", "AB", "SK", "ON", "QC"))+#,
  geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
  scale_y_continuous(name="",expand = c(0,0),
                     limits = c(0,75),
                     breaks=c(0,25,50,75), labels=c("0%", "25%","50%", "75%")) + # controls the y axis
  scale_x_continuous(name="",breaks=c(1,3,6,10,14), labels=c("W1 (Feb)", "W2 (Apr)", "W3 (Jul)", "W4 (Nov)", "W5 (Apr)")) + # controls the x axis, including year ticks
  geom_vline(xintercept = 1.75, linetype="dotted")+
  geom_vline(xintercept = 4, linetype="solid")+
  geom_vline(xintercept = 9.5, linetype="dashed")+
  theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
        legend.text = element_text(size=9), # size of legend text
        legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
        legend.key=element_blank(), 
        aspect.ratio = 0.75, # of the entire grid (y vs. x axes)
        axis.text=element_text(size=9), # of the year ticks on the x axis
        axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
  
  theme( # remove the vertical grid lines
    panel.grid=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"))+ # add x and y axis lines
  geom_hline(yintercept=50,colour="gray")

dev.copy2pdf(file="PaperFiguresFinal/strongoppose_overtime5.pdf")



##### SI SECTION 5: DIFFERENTIAL TRENDS IN CARBON PRICING SUPPORT ACROSS REBATE AND NON-REBATE PROVINCES ##### 

# For parallel trends figures (Figure A3 and Figure A4)
pollingdata_canada <- read.csv("ReplicationData_CanadianSurveysonEnergyandEnvironment.csv")

####### 

# first for carbont tax question
taxsupport_summary <- pollingdata_canada %>%
  group_by(year) %>%
  dplyr::summarize(taxBC = mean(carbon_tax_bin[which(provcode=="BC")], na.rm=TRUE),
                   taxAB = mean(carbon_tax_bin[which(provcode=="AB")], na.rm=TRUE),
                   taxON = mean(carbon_tax_bin[which(provcode=="ON")], na.rm=TRUE),
                   taxQC = mean(carbon_tax_bin[which(provcode=="QC")], na.rm=TRUE),
                   seBC = sqrt(var(carbon_tax_bin[which(provcode=="BC")], na.rm=TRUE)/length(carbon_tax_bin[which(provcode=="BC")])),
                   seAB = sqrt(var(carbon_tax_bin[which(provcode=="AB")], na.rm=TRUE)/length(carbon_tax_bin[which(provcode=="AB")])),
                   seON = sqrt(var(carbon_tax_bin[which(provcode=="ON")], na.rm=TRUE)/length(carbon_tax_bin[which(provcode=="ON")])),
                   seQC = sqrt(var(carbon_tax_bin[which(provcode=="QC")], na.rm=TRUE)/length(carbon_tax_bin[which(provcode=="QC")])))

tax_supportmelt <- melt(taxsupport_summary, id="year")
tax_supportmelt$year<- as.numeric(as.character(tax_supportmelt$year))
tax_supportmelt$se <- tax_supportmelt$value[21:40]

### Question wording is different across waves for carbon tax
#Consistent Wording A: 2011, 2013, 2014
# 2015 is inconsistent, as is 2011 for Quebec
# 2018 is not quite the same but sufficiently close 

# Drop inconsistent wording
tax_supportmelt <- tax_supportmelt[-which(tax_supportmelt$year==2015),] 
tax_supportmelt <- tax_supportmelt[-which(tax_supportmelt$year==2011&tax_supportmelt$variable=="taxQC"),] # also different
tax_supportmelt <- tax_supportmelt[1:15,]

quartz(width=7, height=5.5)

ggplot(tax_supportmelt, aes(x = year, y = value*100)) + 
  ggtitle("Support for carbon taxation, by year and province") + # title
  geom_line(aes(color = variable), size = linesize) + # aesthetics of the line between points
  geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
  scale_color_manual(name="",
                     values=c("#E69F00", "#56B4E9", "#D55E00","#CC79A7"),#, , "#D55E00", "#CC79A7")
                     labels=c("BC", "AB","ON", "QC"))+#,
  geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
  scale_y_continuous(name="",expand = c(0,0),
                     limits = c(0,100),
                     breaks=c(0,25,50,75,100), labels=c("0%", "25%","50%", "75%", "100%")) + # controls the y axis
  scale_x_continuous(name="")+#,breaks=c(1,3,6), labels=c("Wave 1", "Wave 2", "Wave3")) + # controls the x axis, including year ticks
  theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
        legend.text = element_text(size=9), # size of legend text
        legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
        legend.key=element_blank(), 
        axis.text=element_text(size=9), # of the year ticks on the x axis
        axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
  
  theme( # remove the vertical grid lines
    panel.grid=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"))+ # add x and y axis lines
  geom_hline(yintercept=50,colour="gray")

dev.copy2pdf(file="PaperFiguresFinal/tax_overtime_wordingA.pdf")


# second, for cap and trade question
capsupport_summary <- pollingdata_canada %>%
  group_by(year) %>%
  dplyr::summarize(capBC = mean(captrade_bin[which(provcode=="BC")], na.rm=TRUE),
                   capAB = mean(captrade_bin[which(provcode=="AB")], na.rm=TRUE),
                   capON = mean(captrade_bin[which(provcode=="ON")], na.rm=TRUE),
                   capQC = mean(captrade_bin[which(provcode=="QC")], na.rm=TRUE),
                   seBC = sqrt(var(captrade_bin[which(provcode=="BC")], na.rm=TRUE)/length(captrade_bin[which(provcode=="BC")])),
                   seAB = sqrt(var(captrade_bin[which(provcode=="AB")], na.rm=TRUE)/length(captrade_bin[which(provcode=="AB")])),
                   seON = sqrt(var(captrade_bin[which(provcode=="ON")], na.rm=TRUE)/length(captrade_bin[which(provcode=="ON")])),
                   seQC = sqrt(var(captrade_bin[which(provcode=="QC")], na.rm=TRUE)/length(captrade_bin[which(provcode=="QC")])))

cap_supportmelt <- melt(capsupport_summary, id="year")
cap_supportmelt$year<- as.numeric(as.character(cap_supportmelt$year))
cap_supportmelt$se <- cap_supportmelt$value[21:40]

cap_supportmelt<- cap_supportmelt[-which(cap_supportmelt$year==2013&cap_supportmelt$variable=="capQC"),]
cap_supportmelt <- cap_supportmelt[1:19,]


quartz(width=7, height=5.5)

ggplot(cap_supportmelt, aes(x = year, y = value*100)) + 
  ggtitle("Support for emissions trading, by province and year") + # title
  geom_line(aes(color = variable), size = linesize) + # aesthetics of the line between points
  geom_point(aes(color = variable), size = pointsize)+ # aesthetics of the points
  scale_color_manual(name="",
                     values=c("#E69F00", "#56B4E9", "#D55E00","#CC79A7"),#, , "#D55E00", "#CC79A7")
                     labels=c("BC", "AB","ON", "QC"))+#,
  geom_errorbar(aes(ymin=value*100-se*100, ymax=value*100+se*100,color=variable), width=.2)+
  scale_y_continuous(name="",expand = c(0,0),
                     limits = c(0,100),
                     breaks=c(0,25,50,75,100), labels=c("0%", "25%","50%", "75%", "100%")) + # controls the y axis
  scale_x_continuous(name="")+#,breaks=c(1,3,6), labels=c("Wave 1", "Wave 2", "Wave3")) + # controls the x axis, including year ticks
  theme(plot.title = element_text(size = 12, face = "bold"),#, hjust=0.5), # controls title font + size
        legend.text = element_text(size=9), # size of legend text
        legend.key.size = unit(1.1, 'lines'), # shape of legend boxes
        legend.key=element_blank(), 
        axis.text=element_text(size=9), # of the year ticks on the x axis
        axis.text.x=element_text(vjust = 0))+ # controls the orientation of these x axis labels
  
  theme( # remove the vertical grid lines
    panel.grid=element_blank(),
    panel.background = element_blank(),
    panel.border = element_blank(), #element_rect(colour = "black", fill=NA, size=0.25),
    axis.line=element_line(colour="gray"))+ # add x and y axis lines
  geom_hline(yintercept=50,colour="gray")

dev.copy2pdf(file="PaperFiguresFinal/cap_overtime.pdf")



# Table A1
# Wave 2 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta, all respondents
did.reg=lm(support_cp~postperiod_1_2*fedprice, data = cpanel[-which(cpanel$wave=="wave3"|cpanel$wave=="wave4"|cpanel$wave=="wave5"|cpanel$prov_clean=="AB"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 2 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta, only repeat respondents through Wave 2
did.reg=lm(support_cp~postperiod_1_2*fedprice, data = cpanel.complete2[-which(cpanel.complete2$wave=="wave3"|
                                                                                cpanel.complete2$wave=="wave4"|cpanel.complete2$wave=="wave5"|cpanel.complete2$prov_clean=="AB"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 2 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta, only repeat respondents through Wave 5
did.reg=lm(support_cp~postperiod_1_2*fedprice, data = cpanel.complete5[-which(cpanel.complete5$wave=="wave3"|
                                                                                cpanel.complete5$wave=="wave4"|cpanel.complete5$wave=="wave5"|cpanel.complete5$prov_clean=="AB"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))


# Wave 3 vs Wave 2. DV: Liklihood of supporting CP policy, dropping Alberta, all respondents
did.reg=lm(support_cp~postperiod_2_3*fedprice, data = cpanel[-which(cpanel$wave=="wave1"|cpanel$wave=="wave4"|cpanel$wave=="wave5"|cpanel$prov_clean=="AB"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 3 vs Wave 2. DV: Liklihood of supporting CP policy, dropping Alberta, only repeat respondents through Wave 3
did.reg=lm(support_cp~postperiod_2_3*fedprice, data = cpanel.complete2[-which(cpanel.complete3$wave=="wave1"|
                                                                                cpanel.complete3$wave=="wave4"|cpanel.complete3$wave=="wave5"|cpanel.complete3$prov_clean=="AB"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 3 vs Wave 2. DV: Liklihood of supporting CP policy, dropping Alberta, only repeat respondents through Wave 5
did.reg=lm(support_cp~postperiod_2_3*fedprice, data = cpanel.complete5[-which(cpanel.complete5$wave=="wave1"|
                                                                                cpanel.complete5$wave=="wave4"|cpanel.complete5$wave=="wave5"|cpanel.complete5$prov_clean=="AB"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))


# Wave 5 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta, all respondents
did.reg=lm(support_cp~postperiod_1_5*fedprice, data = cpanel[-which(cpanel$wave=="wave2"|cpanel$wave=="wave3"|cpanel$wave=="wave4"|cpanel$prov_clean=="AB"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 5 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta, only repeat respondents through Wave 5
did.reg=lm(support_cp~postperiod_1_5*fedprice, data = cpanel.complete5[-which(cpanel.complete5$wave=="wave2"|
                                                                                cpanel.complete5$wave=="wave3"|cpanel.complete5$wave=="wave4"|cpanel.complete5$prov_clean=="AB"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))



## Table A2 (dropping Ontario)

# Wave 2 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta + Ontario, all respondents
did.reg=lm(support_cp~postperiod_1_2*fedprice, data = cpanel[-which(cpanel$wave=="wave3"|cpanel$wave=="wave4"|cpanel$wave=="wave5"|cpanel$prov_clean=="AB"|cpanel$prov_clean=="ON"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 2 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta + Ontario, only repeat respondents through Wave 2
did.reg=lm(support_cp~postperiod_1_2*fedprice, data = cpanel.complete2[-which(cpanel.complete2$wave=="wave3"|
                                                                                cpanel.complete2$wave=="wave4"|cpanel.complete2$wave=="wave5"|cpanel.complete2$prov_clean=="AB"|cpanel.complete2$prov_clean=="ON"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 2 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta + Ontario, only repeat respondents through Wave 5
did.reg=lm(support_cp~postperiod_1_2*fedprice, data = cpanel.complete5[-which(cpanel.complete5$wave=="wave3"|
                                                                                cpanel.complete5$wave=="wave4"|cpanel.complete5$wave=="wave5"|cpanel.complete5$prov_clean=="AB"|cpanel.complete5$prov_clean=="ON"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))


# Wave 3 vs Wave 2. DV: Liklihood of supporting CP policy, dropping Alberta + Ontario, all respondents
did.reg=lm(support_cp~postperiod_2_3*fedprice, data = cpanel[-which(cpanel$wave=="wave1"|cpanel$wave=="wave4"|cpanel$wave=="wave5"|cpanel$prov_clean=="AB"|cpanel$prov_clean=="ON"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 3 vs Wave 2. DV: Liklihood of supporting CP policy, dropping Alberta + Ontario, only repeat respondents through Wave 3
did.reg=lm(support_cp~postperiod_2_3*fedprice, data = cpanel.complete2[-which(cpanel.complete3$wave=="wave1"|
                                                                                cpanel.complete3$wave=="wave4"|cpanel.complete3$wave=="wave5"|cpanel.complete3$prov_clean=="AB"|cpanel.complete3$prov_clean=="ON"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 3 vs Wave 2. DV: Liklihood of supporting CP policy, dropping Alberta + Ontario, only repeat respondents through Wave 5
did.reg=lm(support_cp~postperiod_2_3*fedprice, data = cpanel.complete5[-which(cpanel.complete5$wave=="wave1"|
                                                                                cpanel.complete5$wave=="wave4"|cpanel.complete5$wave=="wave5"|cpanel.complete5$prov_clean=="AB"|cpanel.complete5$prov_clean=="ON"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))


# Wave 5 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta + Ontario, all respondents
did.reg=lm(support_cp~postperiod_1_5*fedprice, data = cpanel[-which(cpanel$wave=="wave2"|cpanel$wave=="wave3"|cpanel$wave=="wave4"|cpanel$prov_clean=="AB"|cpanel$prov_clean=="ON"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))

# Wave 5 vs Wave 1. DV: Liklihood of supporting CP policy, dropping Alberta + Ontario, only repeat respondents through Wave 5
did.reg=lm(support_cp~postperiod_1_5*fedprice, data = cpanel.complete5[-which(cpanel.complete5$wave=="wave2"|
                                                                                cpanel.complete5$wave=="wave3"|cpanel.complete5$wave=="wave4"|cpanel.complete5$prov_clean=="AB"|cpanel.complete5$prov_clean=="ON"),])
coeftest(did.reg,vcov=function(x)
  vcovHC(x, cluster="responseid", type="HC1"))



##### SI SECTION 6: CARBON PRICING SUPPORT THROUGH TIME BY PARTY PREFERENCE ##### 

## Table A3
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Liberal Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Conservative Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="ndp")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Green Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="I would not vote")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Liberal Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Conservative Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="ndp")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Green Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave1"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="I would not vote")])

mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Liberal Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Conservative Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="ndp")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Green Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="I would not vote")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Liberal Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Conservative Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="ndp")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Green Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave2"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="I would not vote")])

mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Liberal Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Conservative Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="ndp")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Green Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="I would not vote")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Liberal Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Conservative Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="ndp")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Green Party")])
mean(cpanel$support_cp[which(cpanel$wave=="wave3"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="I would not vote")])

mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Liberal Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Conservative Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="ndp")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Green Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="I would not vote")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Liberal Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Conservative Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="ndp")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Green Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave4"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="I would not vote")], na.rm=TRUE)

mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Liberal Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Conservative Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="ndp")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="Green Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="ON"&cpanel$wave1_party=="I would not vote")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Liberal Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Conservative Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="ndp")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="Green Party")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave5"&cpanel$prov_clean=="SK"&cpanel$wave1_party=="I would not vote")], na.rm=TRUE)



##### SI SECTION 7: No code ##### 

##### SI SECTION 8: No code ##### 


##### SI SECTION 9: PERSISTENCE OF REBATER MISPERCEPTIONS IN WAVE 4 ##### 

prop.table(table(cpanel$think_rebate, cpanel$prov_clean), margin=2) 

mean(cpanel$divestimate[which(cpanel$prov_clean=="BC"&cpanel$wave=="wave4")], na.rm=TRUE)
mean(cpanel$divestimate[which(cpanel$prov_clean=="AB"&cpanel$wave=="wave4")], na.rm=TRUE)
mean(cpanel$divestimate[which(cpanel$prov_clean=="QC"&cpanel$wave=="wave4")], na.rm=TRUE)

sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="BC"&cpanel$wave=="wave4")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="BC"&cpanel$wave=="wave4")]))
sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="AB"&cpanel$wave=="wave4")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="AB"&cpanel$wave=="wave4")]))
sqrt(var(cpanel$divestimate[which(cpanel$prov_clean=="QC"&cpanel$wave=="wave4")], na.rm=TRUE)/length(cpanel$divestimate[which(cpanel$prov_clean=="QC"&cpanel$wave=="wave4")]))



#### SI SECTION 10: DISTRIBUTION OF REBATE PERCEPTIONS AMONG ONTARIO AND SASKATCHEWAN RESIDENTS #### 

# Figure A5
ggplot(divcheck, aes(x=federal_realdividend_raw)) + 
  geom_bar(show.legend=FALSE)+
  scale_y_continuous(limits=c(0,200)) +
  ylab("Number of survey respondents") +
  xlab("Actual household dividend ($) received by panelist between Wave 2 and Wave 3") +
  theme_bw()


# Descriptive statistics reported in text
subset_wave3 <-cpanel[which(cpanel$wave=="wave3"),]
table(subset_wave3$federal_realdividend_raw[which(subset_wave3$prov_clean=="ON")]>subset_wave3$divestimate_low[which(subset_wave3$prov_clean=="ON")]&
        subset_wave3$federal_realdividend_raw[which(subset_wave3$prov_clean=="ON")]<subset_wave3$divestimate_high[which(subset_wave3$prov_clean=="ON")])
73/(73+228)

table(subset_wave3$federal_realdividend_raw[which(subset_wave3$prov_clean=="SK")]>subset_wave3$divestimate_low[which(subset_wave3$prov_clean=="SK")]&
        subset_wave3$federal_realdividend_raw[which(subset_wave3$prov_clean=="SK")]<subset_wave3$divestimate_high[which(subset_wave3$prov_clean=="SK")])
57/(47+247)

table(subset_wave3$federal_realdividend_raw[which(subset_wave3$prov_clean=="ON")]>subset_wave3$divestimate_high[which(subset_wave3$prov_clean=="ON")])
185/(186+116)

table(subset_wave3$federal_realdividend_raw[which(subset_wave3$prov_clean=="SK")]>subset_wave3$divestimate_high[which(subset_wave3$prov_clean=="SK")])
229/(229+75)

table(subset_wave3$federal_realdividend_raw[which(subset_wave3$prov_clean=="ON")]<subset_wave3$divestimate_low[which(subset_wave3$prov_clean=="ON")])
43/(43+258)

table(subset_wave3$federal_realdividend_raw[which(subset_wave3$prov_clean=="SK")]<subset_wave3$divestimate_low[which(subset_wave3$prov_clean=="SK")])
18/(18+286)



# estimates among partisans
mean(subset_wave3$divestimate[which(subset_wave3$prov_clean=="SK"&subset_wave3$votertype2=="Consistent Conservatives")], na.rm=TRUE)
mean(subset_wave3$divestimate[which(subset_wave3$prov_clean=="ON"&subset_wave3$votertype2=="Consistent Conservatives")], na.rm=TRUE)

mean(subset_wave3$divestimate[which(subset_wave3$prov_clean=="SK"&subset_wave3$votertype2=="Consistent Liberals")], na.rm=TRUE)
mean(subset_wave3$divestimate[which(subset_wave3$prov_clean=="ON"&subset_wave3$votertype2=="Consistent Liberals")], na.rm=TRUE)

sqrt(var(subset_wave3$divestimate[which(subset_wave3$prov_clean=="SK"&subset_wave3$votertype2=="Consistent Conservatives")], na.rm=TRUE)/length(subset_wave3$divestimate[which(subset_wave3$prov_clean=="SK"&subset_wave3$votertype2=="Consistent Conservatives")]))
sqrt(var(subset_wave3$divestimate[which(subset_wave3$prov_clean=="ON"&subset_wave3$votertype2=="Consistent Conservatives")], na.rm=TRUE)/length(subset_wave3$divestimate[which(subset_wave3$prov_clean=="ON"&subset_wave3$votertype2=="Consistent Conservatives")]))
sqrt(var(subset_wave3$divestimate[which(subset_wave3$prov_clean=="SK"&subset_wave3$votertype2=="Consistent Liberals")], na.rm=TRUE)/length(subset_wave3$divestimate[which(subset_wave3$prov_clean=="SK"&subset_wave3$votertype2=="Consistent Liberals")]))
sqrt(var(subset_wave3$divestimate[which(subset_wave3$prov_clean=="ON"&subset_wave3$votertype2=="Consistent Liberals")], na.rm=TRUE)/length(subset_wave3$divestimate[which(subset_wave3$prov_clean=="ON"&subset_wave3$votertype2=="Consistent Liberals")]))


t.test(subset_wave3$divestimate[which(subset_wave3$votertype14=="Consistent Conservatives"&subset_wave3$prov_clean=="SK")],
       subset_wave3$divestimate[which(subset_wave3$votertype14=="Consistent Liberals"&subset_wave3$prov_clean=="SK")],
       var.equal=FALSE, na.rm=TRUE)

t.test(subset_wave3$divestimate[which(subset_wave3$votertype14=="Consistent Conservatives"&subset_wave3$prov_clean=="ON")],
       subset_wave3$divestimate[which(subset_wave3$votertype14=="Consistent Liberals"&subset_wave3$prov_clean=="ON")],
       var.equal=FALSE, na.rm=TRUE)





#### SI SECTION 11: DO CANADIANS CORRECTLY REPORT PROVINCIAL VS. FEDERAL TAX REGIME? #### 


#fed price places
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave1"&(cpanel$prov_clean=="ON"|cpanel$prov_clean=="SK"))]))
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave2"&(cpanel$prov_clean=="ON"|cpanel$prov_clean=="SK"))]))
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave3"&(cpanel$prov_clean=="ON"|cpanel$prov_clean=="SK"))]))
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave4"&(cpanel$prov_clean=="ON"|cpanel$prov_clean=="SK"))]))

# prov price places
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave1"&(cpanel$prov_clean=="BC"|cpanel$prov_clean=="QC"))]))
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave2"&(cpanel$prov_clean=="BC"|cpanel$prov_clean=="QC"))]))
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave3"&(cpanel$prov_clean=="BC"|cpanel$prov_clean=="QC"))]))
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave4"&(cpanel$prov_clean=="BC"|cpanel$prov_clean=="QC"))]))

#special case of Alberta
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave1"&cpanel$prov_clean=="AB")]))
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave2"&cpanel$prov_clean=="AB")]))
prop.table(table(cpanel$prov_or_fed[which(cpanel$wave=="wave3"&cpanel$prov_clean=="AB")]))
prop.table(table(cpanel$po7b[which(cpanel$wave=="wave4"&cpanel$prov_clean=="AB")]))


#### SI SECTION 12: No code #### 


#### SI SECTION 13: EXPERIMENTAL BALANCE IN CANADA #### 

t.test(cpanel$genderBIN[which(cpanel$treat=="Control")],
       cpanel$genderBIN[which(cpanel$treat=="Treat")],
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$age[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$age[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$bachelorsBIN[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$bachelorsBIN[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$ideology[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$ideology[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$trust.CRA[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$trust.CRA[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$trust.FED[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$trust.FED[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$own[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$own[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$numbercars[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$numbercars[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$heat_perceived[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$heat_perceived[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)

t.test(as.numeric(cpanel$gas_perceived[which(cpanel$treat=="Control")]),
       as.numeric(cpanel$gas_perceived[which(cpanel$treat=="Treat")]),
       var.equal=FALSE, na.rm=TRUE)



#### SI SECTION 14: MANIPULATION CHECK #### 

treatedgroup <- cpanel[which(cpanel$treat=="Treat"),]
newrespondent_control <- wave4.pool[which(wave4.pool$prov_clean=="ON"|wave4.pool$prov_clean=="SK"),] 

t.test(as.numeric(treatedgroup$divestimate_posttreatment),
       as.numeric(newrespondent_control$perceived_rebatecategory),
       var.equal=FALSE, na.rm=TRUE)

sqrt(var(as.numeric(treatedgroup$divestimate_posttreatment)))
sqrt(var(as.numeric(newrespondent_control$perceived_rebatecategory)))

t.test(as.numeric(treatedgroup$divestimate_posttreatment),
       as.numeric(treatedgroup$perceivedrebate.category_wave3),
       var.equal=FALSE, na.rm=TRUE)

sqrt(var(as.numeric(treatedgroup$divestimate_posttreatment)))
sqrt(var(as.numeric(treatedgroup$perceivedrebate.category_wave3)))


estimatedifferences <- cpanel$dv0[which(cpanel$treat=="Treat")]-cpanel$federal_realdividend_raw[which(cpanel$treat=="Treat")]

length(estimatedifferences)
sum(estimatedifferences< -50)
sum(estimatedifferences> 50)
sum(estimatedifferences< 50 & estimatedifferences > -50)




#### SI SECTION 15: No Code #### 

#### SI SECTION 16: See Swiss replication script #### 

#### SI SECTION 17: No Code #### 

#### SI SECTION 18: PARTISAN DIFFERENCES IN POLITIAL AD RECOLLECTION  #### 


mean(cpanel$adexposure[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")], na.rm=TRUE)
sqrt(var(cpanel$adexposure[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")], na.rm=TRUE))

mean(cpanel$adexposure[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")], na.rm=TRUE)
sqrt(var(cpanel$adexposure[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")], na.rm=TRUE))

t.test(cpanel$adexposure[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")],
       cpanel$adexposure[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")])

mean(cpanel$adnegative[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")], na.rm=TRUE)
sqrt(var(cpanel$adnegative[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")], na.rm=TRUE))

mean(cpanel$adnegative[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")], na.rm=TRUE)
sqrt(var(cpanel$adnegative[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")], na.rm=TRUE))

t.test(cpanel$adnegative[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")],
       cpanel$adnegative[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")])

mean(cpanel$adinfluence[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")], na.rm=TRUE)
sqrt(var(cpanel$adinfluence[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")], na.rm=TRUE))

mean(cpanel$adinfluence[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")], na.rm=TRUE)
sqrt(var(cpanel$adinfluence[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")], na.rm=TRUE))

t.test(cpanel$adinfluence[which(cpanel$votertype34=="Consistent Conservatives"&cpanel$wave=="wave4")],
       cpanel$adinfluence[which(cpanel$votertype34=="Consistent Liberals"&cpanel$wave=="wave4")])


#### SI SECTION 19: PARTISAN ACCURACY OF POST-TREATMENT REBATE ESTIMATES #### 

exp_treatedgroup <- cpanel[which(cpanel$treat=="Treat"),] 

exp_treatedgroup$estimate_difference <- exp_treatedgroup$federal_realdividend_raw - exp_treatedgroup$dv0
min(exp_treatedgroup$estimate_difference)
max(exp_treatedgroup$estimate_difference)

mean(exp_treatedgroup$estimate_difference[which(exp_treatedgroup$votertype34=="Consistent Conservatives")])
sqrt(var(exp_treatedgroup$estimate_difference[which(exp_treatedgroup$votertype34=="Consistent Conservatives")]))

mean(exp_treatedgroup$estimate_difference[which(exp_treatedgroup$votertype34=="Consistent Liberals")])
sqrt(var(exp_treatedgroup$estimate_difference[which(exp_treatedgroup$votertype34=="Consistent Liberals")]))


t.test(exp_treatedgroup$estimate_difference[which(exp_treatedgroup$votertype34=="Consistent Conservatives")],
       exp_treatedgroup$estimate_difference[which(exp_treatedgroup$votertype34=="Consistent Liberals")])






#### SI SECTION 20: EXPLORING SURVEY DESIGN EFFECTS #### 
mean(cpanel$heardBIN_atall[which(cpanel$wave=="wave2")], na.rm=TRUE)
mean(cpanel$heardBIN_atall[which(cpanel$wave=="wave3")], na.rm=TRUE)
mean(cpanel$heardBIN_atall[which(cpanel$wave=="wave4")], na.rm=TRUE)

mean(wave2.pool$heardBIN_atall, na.rm=TRUE)
mean(wave3.pool$heardBIN_atall, na.rm=TRUE)
mean(wave4.pool$heardBIN_atall, na.rm=TRUE)


mean(cpanel$heardBIN_modormore[which(cpanel$wave=="wave2")], na.rm=TRUE)
mean(cpanel$heardBIN_modormore[which(cpanel$wave=="wave3")], na.rm=TRUE)
mean(cpanel$heardBIN_modormore[which(cpanel$wave=="wave4")], na.rm=TRUE)

mean(wave2.pool$heardBIN_modormore, na.rm=TRUE)
mean(wave3.pool$heardBIN_modormore, na.rm=TRUE)
mean(wave4.pool$heardBIN_modormore, na.rm=TRUE)


mean(cpanel$support_cp[which(cpanel$wave=="wave2")], na.rm=TRUE)
mean(cpanel$support_cp[which(cpanel$wave=="wave3")], na.rm=TRUE)

mean(wave2.pool$support_cp, na.rm=TRUE)
mean(wave3.pool$support_cp, na.rm=TRUE)
mean(wave4.pool$support_cp, na.rm=TRUE)


mean(cpanel$fair_cp[which(cpanel$wave=="wave2")], na.rm=TRUE)
mean(cpanel$fair_cp[which(cpanel$wave=="wave3")], na.rm=TRUE)
mean(cpanel$fair_cp[which(cpanel$wave=="wave4")], na.rm=TRUE)

mean(wave2.pool$fair_cp, na.rm=TRUE)
mean(wave3.pool$fair_cp, na.rm=TRUE)
mean(wave4.pool$fair_cp, na.rm=TRUE)


mean(cpanel$correct_p_or_f[which(cpanel$wave=="wave2")], na.rm=TRUE)
mean(cpanel$correct_p_or_f[which(cpanel$wave=="wave3")], na.rm=TRUE)
mean(cpanel$correct_p_or_f[which(cpanel$wave=="wave4")], na.rm=TRUE)

mean(wave2.pool$correct_p_or_f, na.rm=TRUE)
mean(wave3.pool$correct_p_or_f, na.rm=TRUE)
mean(wave4.pool$correct_p_or_f, na.rm=TRUE)


#### SI SECTION 21: SAMPLE REPRESENTATIVENESS IN CANADA #### 

wave1sub <- cpanel[which(cpanel$wave=="wave1"),]

round(prop.table(table(wave1sub$gender, wave1sub$prov_clean), margin=2), digits=4) 
round(prop.table(table(wave1sub$frenchBIN, wave1sub$prov_clean), margin=2), digits=4) 
round(prop.table(table(wave1sub$age3, wave1sub$prov_clean), margin=2), digits=4) 
round(prop.table(table(wave1sub$edu5, wave1sub$prov_clean), margin=2), digits=4) 
round(prop.table(table(wave1sub$income6, wave1sub$prov_clean), margin=2), digits=4) 



#### SI SECTION 22: PANEL ATTRITION #### 

attrit12 <- wave1_IDs[!(wave1_IDs%in%wave2_IDs)]
attrit23 <- wave2_IDs[!(wave2_IDs%in%wave3_IDs)]
attrit34 <- wave3_IDs[!(wave3_IDs%in%wave4_IDs)]
attrit45 <- wave4_IDs[!(wave4_IDs%in%wave5_IDs)]

# Table A11
length(wave1_IDs)
length(attrit12)
length(wave2_IDs)
length(attrit23)
length(wave3_IDs)
length(attrit34)
length(wave4_IDs)
length(attrit45)
length(wave5_IDs)


# Table A12
#gender
mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave1_IDs], na.rm=TRUE)
mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit12], na.rm=TRUE)

t.test(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave1_IDs],
       cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit12],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave2_IDs], na.rm=TRUE)
mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit23], na.rm=TRUE)

t.test(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave2_IDs],
       cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit23],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave3_IDs], na.rm=TRUE)
mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit34], na.rm=TRUE)

t.test(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave3_IDs],
       cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit34],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave4_IDs], na.rm=TRUE)
mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit45], na.rm=TRUE)

t.test(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave4_IDs],
       cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit45],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$genderBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave5_IDs], na.rm=TRUE)


#language
mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave1_IDs], na.rm=TRUE)
mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit12], na.rm=TRUE)

t.test(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave1_IDs],
       cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit12],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave2_IDs], na.rm=TRUE)
mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit23], na.rm=TRUE)

t.test(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave2_IDs],
       cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit23],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave3_IDs], na.rm=TRUE)
mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit34], na.rm=TRUE)

t.test(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave3_IDs],
       cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit34],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave4_IDs], na.rm=TRUE)
mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit45], na.rm=TRUE)

t.test(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave4_IDs],
       cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit45],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$frenchBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave5_IDs], na.rm=TRUE)


#income
mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave1_IDs], na.rm=TRUE)
mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%attrit12], na.rm=TRUE)

t.test(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave1_IDs],
       cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%attrit12],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave2_IDs], na.rm=TRUE)
mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%attrit23], na.rm=TRUE)

t.test(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave2_IDs],
       cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%attrit23],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave3_IDs], na.rm=TRUE)
mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%attrit34], na.rm=TRUE)

t.test(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave3_IDs],
       cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%attrit34],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave4_IDs], na.rm=TRUE)
mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%attrit45], na.rm=TRUE)

t.test(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave4_IDs],
       cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%attrit45],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$income.NUM[cpanel$wave=="wave1"&cpanel$responseid%in%wave5_IDs], na.rm=TRUE)


#income
mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave1_IDs], na.rm=TRUE)
mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit12], na.rm=TRUE)

t.test(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave1_IDs],
       cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit12],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave2_IDs], na.rm=TRUE)
mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit23], na.rm=TRUE)

t.test(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave2_IDs],
       cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit23],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave3_IDs], na.rm=TRUE)
mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit34], na.rm=TRUE)

t.test(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave3_IDs],
       cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit34],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave4_IDs], na.rm=TRUE)
mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit45], na.rm=TRUE)

t.test(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave4_IDs],
       cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%attrit45],
       var.equal=FALSE, na.rm=TRUE)

mean(cpanel$bachelorsBIN[cpanel$wave=="wave1"&cpanel$responseid%in%wave5_IDs], na.rm=TRUE)


#### SI SECTION 23: No Code #### 

#### SI SECTION 24: No Code #### 

#### SI SECTION 25: See Swiss Replication Script #### 

#### SI SECTION 26: No Code #### 

#### SI SECTION 27: See Swiss Replication Script #### 















